# Bibliographic Information & Notes ############################################
# Title: Taking the Cloth -- Social norms and elite cues increase support for... 
# ...masks among white Evangelical Americans
# By: Claire Adida, Christina Cottiero, Leonardo Falabella, Isabel Gotti, ...
# ...ShahBano Ijaz, Gregoire Phillips, and Michael Seese
# Published In: Journal of Experimental Political Science
# Publication Date: 
# DOI: 
# Preprint: 10.31235/osf.io/yt3e7
# Note: This code replicates findings in the published manuscript, for SI...
# ...materials, please see accompanying file TtC_SI_Replication.R


# Packages and Data ############################################################
## Working Directory ###########################################################
rm(list = ls())
setwd("/Users/mikeseese/Documents/GitHub/Taking_the_Cloth_Replication") # Change WD


## Packages ####################################################################
library("tidyverse")
library("modelsummary")
library("ggeffects")
library("esc")
library("factoextra")


## Data ########################################################################
working <- readRDS(file = "working.rds")


# Preliminary Data Cleaning ####################################################
## Recode Religiosity Index ####################################################
## Corrects a data transfer error from Lucid
working$rel1x <- dplyr::recode(working$rel1,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel2x <- dplyr::recode(working$rel2,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel3x <- dplyr::recode(working$rel3,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel.ix <- working$rel1x + working$rel2x + working$rel3x


# Descriptive Statistics #######################################################
## Descriptives: Covariates ####################################################
x1 <- working %>% select(
  `Age` = age,
  `Follows News` = news.follow,
  `Follows Covid` = covid.follow,
  `Covid Threat` = covid.threat)

datasummary_skim(x1)

x2 <- working %>% select(
  `Gender` = gender,
  `Education` = edu,
  `Party ID` = party.id,
  `News Source` = news.source,
  `Covid News Source` = covid.source,
  `Tested for Covid` = covid.tested,
  `Friends with Covid` = covid.friends)

datasummary_skim(x2, type = "categorical")


## Descriptives: Religious Variables ###########################################
x3 <- working %>% select(
  `Regularly Attends Church` = rel1x,
  `Spiritual Values` = rel2x,
  `Americans Should be More Religious` = rel3x,
  `Religiosity Index` = rel.ix,
  `Religiosity (Self-Reported)` = religiosity,
  `Religious Motivation - Insurance` = c1, 
  `Religious Motivation - Morality` = c2, 
  `Religious Motivation - Social Capital` = c3, 
  `Religious Motivation - Heuristics` = c4,
  `Religious Motivation - Tradition` = c5, 
  `Religious Motivation - Social Identity` = c6, 
  `Religion Offers Comfort` = ei1,
  `Takes Religious Approach to Life` = ei2,
  `Enjoys Social Aspects of Church` = ei3,
  `Extrinsic / Intrinsic Motivation Index` = ei.index)

datasummary_skim(x3)

rm(x1, x2, x3, summaryname)


## Gender, Age, and Partisanship ###############################################
## Gender
ggplot(data = working, aes(x = factor(gender))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  labs(title = "Gender Distribution of Sample",
       x = "Gender",
       y = "Frequency") +
  theme_bw()

## Age
ggplot(data = working, aes(x = age)) + 
  geom_histogram(aes(y = ..density..)) +
  geom_vline(aes(xintercept = mean(age, na.rm = TRUE)),
             color ="red", linetype = "dashed", size = 0.5) +
  labs(title = "Age Distribution of Sample",
       x = "Age",
       y = "Density",
       caption = "Mean Age = 54.23") +
  theme_bw()

## Party ID
ggplot(data = working, aes(x = factor(party.id))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  labs(title = "Partisan Distribution of Sample", 
       subtitle = "Seven Categories",
       x = "Party ID",
       y = "Frequency") +
  theme_bw()

ggplot(data = working, aes(x = factor(party.3))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  labs(title = "Partisan Distribution of Sample", 
       subtitle = "Three Categories",
       x = "Party ID",
       y = "Frequency") +
  theme_bw()

ggplot(data = working, aes(x = factor(party.3))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = 1.5, 
            color = "white") +
  labs(title = "Partisan Distribution of Sample, by Treatment Condition", 
       subtitle = "Three Categories",
       x = "Party ID",
       y = "Frequency") +
  facet_grid(rows = vars(treat)) +
  theme_bw()


# Correlates of Mask Wearing (Control Condition) ###############################
## Demographic Correlates of Mask Wearing ######################################
### Data Cleaning ##############################################################
control <- working %>%
  subset(treat == "Control")

control$ag <- factor(control$age.group, ordered = FALSE)
control$eg <- factor(control$edu, ordered = FALSE)


### OLS Models #################################################################
## OLS Model: Demographics
x1 <- lm(therm ~ gender + ag + eg + party.f + news.follow + covid.follow + 
           covid.threat + covid.tested + covid.friends, data = control)
summary(x1)

## OLS Table
cm <- c("(Intercept)" = "Intercept",
        "genderMale" = "Male",
        "ag30-39" = "Age 30-39",
        "ag40-49" = "Age 40-49",
        "ag50-59" = "Age 50-59",
        "ag60-69" = "Age 60-69",
        "ag70+" = "Age 70+",
        "egHigh School" = "High School",
        "egSome College" = "Some College",
        "egVocational School" = "Vocational School",
        "egAssociate's" = "Associate's",
        "egBachelor's" = "Bachelor's",
        "egSome Graduate School" = "Some Graduate School",
        "egMaster's" = "Master's",
        "egProfessional Degree" = "Professional Degree",
        "egDoctorate" = "Doctorate",
        "party.fDemocrat" = "Democrat",
        "party.fLean Democrat" = "Lean Democrat",
        "party.fIndependent" = "Independent",
        "party.fLean Republican" = "Lean Republican",
        "party.fRepublican" = "Republican",
        "party.fStrong Republican" = "Strong Republican",
        "news.follor" = "Follows News",
        "covid.threat" = "Covid Threat",
        "covid.testedYes" = "Tested for Covid",
        "covid.friendsYes" = "Friends with Covid")

modelsummary(x1,
             # output = "latex",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

rm("cm")

## Linear Predictions by Demographic
ggplot(ggpredict(x1, terms = "ag"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Age Group", y = "Predicted Value")

ggplot(ggpredict(x1, terms = "eg"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Education", y = "Predicted Value")

ggplot(ggpredict(x1, terms = "party.f"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party", y = "Predicted Value")

rm("x1")

## Religious Correlates of Mask Wearing ########################################
### Data Cleaning ##############################################################
## Extrinsic / Intrinsic Motivations
control$ei1.f <- factor(control$ei1)
control$ei2.f <- factor(control$ei2)
control$ei3.f <- factor(control$ei3)

control$ei1.f <- dplyr::recode(control$ei1.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

control$ei2.f <- dplyr::recode(control$ei2.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

control$ei3.f <- dplyr::recode(control$ei3.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree")

## Religious Motivations
control$c.all <- NA
control$c.all[which(control$c1 == 1)] <- 1
control$c.all[which(control$c2 == 1)] <- 2
control$c.all[which(control$c3 == 1)] <- 3
control$c.all[which(control$c4 == 1)] <- 4
control$c.all[which(control$c5 == 1)] <- 5
control$c.all[which(control$c6 == 1)] <- 6

table(control$c.all)
control$c.all <- factor(control$c.all)

## Religiosity
control$religiosity.f <- factor(control$religiosity)
control$religiosity.f <- dplyr::recode(control$religiosity.f,
                                       "1" = "Anti-Religious",
                                       "2" = "Not at all Religious",
                                       "3" = "Slightly Religious",
                                       "4" = "Moderately Religious",
                                       "5" = "Very Religious")

## Religiosity Index
control$rel1.f <- factor(control$rel1)
control$rel2.f <- factor(control$rel2)
control$rel3.f <- factor(control$rel3)

control$rel1.f <- dplyr::recode(control$rel1.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")

control$rel2.f <- dplyr::recode(control$rel2.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")

control$rel3.f <- dplyr::recode(control$rel3.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")


### OLS Models #################################################################
## OLS Model: Extrinsic / Intrinsic Motivations
x1 <- lm(therm ~ ei1.f + ei2.f + ei3.f, data = control)
summary(x1)

## OLS Model: Religious Motivations
x2 <- lm(therm ~ c.all, data = control)
summary(x2)

## OLS Model: Religiosity
x3 <- lm(therm ~ religiosity.f, data = control)
summary(x3)

## OLS Model: Religiosity Index
x4 <- lm(therm ~ rel1.f + rel2.f + rel3.f, data = control)
summary(x4)

## OLS Model: All Religious Variables
x5 <- lm(therm ~ ei1.f + ei2.f + ei3.f + c.all + religiosity.f + rel1.f + 
           rel2.f + rel3.f, data = control)
summary(x5)

## OLS Table
cm <- c("(Intercept)" = "Intercept",
        "ei1.fDisagree" = "Religion Offers Comfort: Disagree",
        "ei1.fNeither Agree nor Disagree" = "Religion Offers Comfort: Neither Agree nor Disagree",
        "ei1.fAgree" = "Religion Offers Comfort: Agree",
        "ei1.fStrongly Agree" = "Religion Offers Comfort: Strongly Agree",
        "ei2.fDisagree" = "Religious Approach to Life: Disagree",
        "ei2.fNeither Agree nor Disagree" = "Religious Approach to Life: Neither Agree nor Disagree",
        "ei2.fAgree" = "Religious Approach to Life: Agree",
        "ei2.fStrongly Agree" = "Religious Approach to Life: Strongly Agree",
        "ei3.fDisagree" = "Enjoys Social Aspects of Church: Disagree",
        "ei3.fNeither Agree nor Disagree" = "Enjoys Social Aspects of Church: Neither Agree nor Disagree",
        "ei3.fAgree" = "Enjoys Social Aspects of Church: Agree",
        "ei3.fStrongly Agree" = "Enjoys Social Aspects of Church: Strongly Agree",
        "c.all2" = "Morality",
        "c.all3" = "Social Capital",
        "c.all4" = "Heuristics",
        "c.all5" = "Tradition",
        "c.all6" = "Social Identity",
        "religiosity.fSlightly Religious" = "Religiosity (SR) Slightly Religious",
        "religiosity.fModerately Religious" = "Religiosity (SR) Moderately Religious",
        "religiosity.fVery Religious" = "Religiosity (SR) Very Religious",
        "rel1.fStrongly Agree" = "Regularly Attends Church: Strongly Agree",
        "rel1.fDisagree" = "Regularly Attends Church: Disagree",
        "rel1.fSomewhat Disagree" = "Regularly Attends Church: Somewhat Disagree",
        "rel1.fSomewhat Agree" = "Regularly Attends Church: Somewhat Agree",
        "rel1.fAgree" = "Regularly Attends Church: Agree",
        "rel2.fStrongly Agree" = "Spiritual Values: Strongly Agree",
        "rel2.fDisagree" = "Spiritual Values: Disagree",
        "rel2.fSomewhat Disagree" = "Spiritual Values: Somewhat Disagree",
        "rel2.fSomewhat Agree" = "Spiritual Values: Somewhat Agree",
        "rel2.fAgree" = "Spiritual Values: Agree",
        "rel3.fStrongly Agree" = "Americans Should be More Religious: Strongly Agree",
        "rel3.fDisagree" = "Americans Should be More Religious: Disagree",
        "rel3.fSomewhat Disagree" = "Americans Should be More Religious: Somewhat Disagree",
        "rel3.fSomewhat Agree" = "Americans Should be More Religious: Somewhat Agree",
        "rel3.fAgree" = "Americans Should be More Religious: Agree")

modelsummary(list(x1, x2, x3, x4, x5),
             # output = "latex",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

rm("x1", "x2", "x3", "x4", "x5", "cm")


## Party ID and Mask Wearing ###################################################
### Data Cleaning ##############################################################
control$party.f2 <- factor(control$party.f, levels = c("Strong Democrat", 
                                                       "Democrat", 
                                                       "Lean Democrat",
                                                       "Independent", 
                                                       "Lean Republican", 
                                                       "Republican", 
                                                       "Strong Republican"))

control$party.32 <- factor(control$party.3, levels = c("Democrat",
                                                       "Independent", 
                                                       "Republican"))


### OLS Models #################################################################
## 7 Category Party ID
x1 <- lm(therm ~ party.f2, data = control)
summary(x1)

m.x1 <- ggpredict(x1, terms = "party.f2")

ggplot(m.x1, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party ID", y = "Predicted Value")

## 3 Category Party ID
x2 <- lm(therm ~ party.32, data = control)
summary(x2)

m.x2 <- ggpredict(x2, terms = "party.32")

ggplot(m.x2, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.075) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party ID", y = "Predicted Value")

rm(x1, m.x1, x2, m.x2, control)


# Main Results: Average Treatment Effects ######################################
## Data Cleaning ###############################################################
w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"


## Baseline OLS Models #########################################################
## Models
a1 <- lm(therm ~ treat, data = working)
a2 <- lm(p1.wash ~ treat, data = working)
a3 <- lm(p2.dist ~ treat, data = working)
a4 <- lm(p3.mask ~ treat, data = working)
a5 <- lm(p4.pray ~ treat, data = working)
a6 <- lm(behave ~ treat, data = working)
a7 <- lm(clicked ~ treat, data = w2)

## Create DF for Coefficient Plot
x1 <- tidy(a1, conf.int = TRUE)
x2 <- tidy(a2, conf.int = TRUE)
x3 <- tidy(a3, conf.int = TRUE)
x4 <- tidy(a4, conf.int = TRUE)
x5 <- tidy(a5, conf.int = TRUE)
x6 <- tidy(a6, conf.int = TRUE)
x7 <- tidy(a7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(a1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(a2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(a3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(a4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(a5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(a6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(a7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept = 0, lty = "11", colour = "grey30") +
  geom_errorbar(aes(cond, ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_errorbar(aes(cond, ymin = low, ymax = high), lwd = 1.15, width = 0) +
  geom_point(aes(cond, estimate)) +
  labs(colour = "", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: Baseline OLS Model",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

## Save DF for Combined Plot
y1 <- y
y1$mdl <- "Baseline"

## Cleaning Up
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)


## Demographic OLS Models ######################################################
## Models
b1 <- lm(as.formula(paste("therm", iv1)), data = working)
b2 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
b3 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
b4 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
b5 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
b6 <- lm(as.formula(paste("behave", iv1)), data = working)
b7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

## Create DF for Coefficient Plot
x1 <- tidy(b1, conf.int = TRUE)
x2 <- tidy(b2, conf.int = TRUE)
x3 <- tidy(b3, conf.int = TRUE)
x4 <- tidy(b4, conf.int = TRUE)
x5 <- tidy(b5, conf.int = TRUE)
x6 <- tidy(b6, conf.int = TRUE)
x7 <- tidy(b7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(b1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(b2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(b3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(b4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(b5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(b6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(b7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept = 0, lty = "11", colour="grey30") +
  geom_errorbar(aes(cond, ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_errorbar(aes(cond, ymin = low, ymax = high), lwd = 1.15, width = 0) +
  geom_point(aes(cond, estimate)) +
  labs(colour = "", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Demographic Controls",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

## Save DF for Combined Plot
y2 <- y
y2$mdl <- "Demographic Controls"

## Cleaning Up
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

## OLS Models All Controls #####################################################
## Models
c1 <- lm(as.formula(paste("therm", iv2)), data = working)
c2 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
c3 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
c4 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
c5 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
c6 <- lm(as.formula(paste("behave", iv2)), data = working)
c7 <- lm(as.formula(paste("clicked", iv2)), data = w2)

## Create DF for Coefficient Plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)
x4 <- tidy(c4, conf.int = TRUE)
x5 <- tidy(c5, conf.int = TRUE)
x6 <- tidy(c6, conf.int = TRUE)
x7 <- tidy(c7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(c4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(c5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(c6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(c7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept and Prayer Outcome
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept = 0, lty = "11", colour = "grey30") +
  geom_errorbar(aes(cond, ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_errorbar(aes(cond, ymin = low, ymax = high), lwd = 1.15, width = 0) +
  geom_point(aes(cond, estimate)) +
  labs(colour = "", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Controls",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

## Save DF for Combined Plot
y3 <- y
y3$mdl <- "All Controls"

## Cleaning Up
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)


## Combined Coefficient Plot ###################################################
## Create a DF
y <- rbind(y1, y2, y3)

y$mdl <- factor(y$mdl, 
                levels = c("Baseline", 
                           "Demographic Controls", 
                           "All Controls"), 
                ordered = TRUE)

## The Plot
ggplot(y) +
  geom_hline(yintercept = 0, lty = "11", colour = "grey30") +
  geom_errorbar(aes(cond, ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_errorbar(aes(cond, ymin = low, ymax = high), lwd = 1.15, width = 0) +
  geom_point(aes(cond, estimate)) +
  labs(colour = "", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out), rows = vars(mdl)) + 
  labs(title = "Coefficient Plot: OLS Models",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())


## Effect Sizes for Thermometer Outcome ########################################
summary(c1)

## Norms
esc_B(b = 0.45198, sdy = 2.995311, grp1n = 597, grp2n = 592, es.type = "g")

## Endorsement
esc_B(b = 0.340044, sdy = 2.995311, grp1n = 597, grp2n = 591, es.type = "g")

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# PCA For Footnote 4 ###########################################################
## Simplify the dataset
x <- working %>% select(
  treat, therm, p1.wash, p2.dist, p3.mask, p4.pray, behave)

x2 <- x %>% select(
  "Thermometer" = therm, 
  "Hand Washing" = p1.wash, 
  "Distancing" = p2.dist, 
  "Mask Wearing" = p3.mask, 
  "Prayer" = p4.pray, 
  "Would Click" = behave)

## Run the PCA
pc1 <- prcomp(na.omit(x2), scale = FALSE)
summary(pc1)


## PCA Plots ###################################################################
fviz_eig(pc1)

fviz_pca_var(pc1, 
                     # col.var = "contrib", 
                     repel = TRUE)


fviz_pca_var(pc1, 
             # col.var = "contrib", 
             repel = TRUE)


